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Abstract. The multiconfigurational time-dependent Hartree-Fock equations are discussed and 
solved for a one-dimensional model of the Helium atom. Results for the ground state energy and 
two-particle density as well as the absorption spectrum are presented and compared to direct 
solutions of the time-dependent Schrodinger equation. 



1. Introduction 

In recent years new powerful radiation sources became available for the precise investigation 
of photoionization processes of matter. New methods have made it possible to observe the 
electronic motion in a time-resolved fashion on a scale of attoseconds [T| [2]. With this, several 
effects, such as strong-field tunneling [3] or time-resolved Auger decay |3j could be studied in 
detail. The explanation of the arising effects is a challenge for theoretical physicists, which 
need to face an old problem, namely the solution of the electronic Schrodinger equation with its 
exponentially growing effort with increasing number of the degrees of freedom. 
Several methods have been developed to circumvent this fundamental limitation. Among 
them are time-dependent density functional theory (TDDFT), see e.g. [5j, the method of 
nonequilibrium Green functions (NEGF), e.g. [SHE], or time-dependent reduced density-matrix 
theory (TDRDM), e.g. [9Hllj. which all aim at projecting the Schrodinger equation on a more 
convenient set of equations requiring only a polynomially growing effort in solution. Despite the 
indisputable successes of these methods, they lack a systematical and practically feasible way 
to achieve convergence to the exact result [j] 

In this paper, we apply a method which provides this mentioned feature, namely time- 
dependent Multiconfigurational Hartree-Fock (MCTDHF). It can either be seen as an extension 
of Hartree-Fock, which includes several Slater determinants (or permanents in the case of 
Bosons) instead of a single one, or as an extension of Configuration Interaction, that employs 
time-dependent single-particle orbitals instead of a fixed basis. In the case of infinitely many 
determinants, the results essentially become exact. In MCTDHF, the exponential problem is 
not really avoided, but it is postponed to much larger systems than in direct solutions of the 
Schrodinger equation. Thus, MCTDHF is applicable to few particle systems of roughly ten 

1 In TDDFT, this is actually the main problem, since the result crucially depends on the choice of the exchange- 
correlation functional. NEGF (and TDRDM) are principally exact, if all equations in a hierarchy of equations 
were taken into account, resp. if all self-energy diagrams were summed up. In practice, however, the hierarchy is 
decoupled already on a low level. 



particles. To become familiar with the method, in this work we consider a directly solvable 
one-dimensional model of the Helium atom, and compare the MCTDHF results with those from 
the time-dependent Schrodinger equation (TDSE). 

The outline of the paper is as follows: After this introduction, we give an overview on 
the MCTDHF formalism, recapitulate the working equations and provide the main ideas of 
our implementation. Subsequently, groundstate as well as time-dependent results for one- 
dimensional Helium are presented. In this paper, we apply the notation of Ref. |12j . 

2. The MCTDHF method 

The system of our interest are few electron atoms in an external electromagnetic field described 
by the Hamiltonian (in atomic units) 

* = E{£-H-«w-'4 + i£h^r (1) 
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where the motion of the nucleus is neglected. In this paper, we focus on the one-dimensional 
Helium atom (Z = 2), for which the singularities in the Coulomb potential are softened by 
a cutoff, see Sec. |3.1| However, since the following theoretical considerations are completely 
general, we use the three-dimensional formulation. The studied systems are assumed to be 
initially in the groundstate, so a wavefunction treatment is appropriate. 

2.1. Overview 

Standard wavepacket propagation methods, e.g. time-dependent Configuration Interaction, 
typically approximate the many-body wavefunction as a linear expansion in a set of basis vectors 
of the subspace %{2M, N) of the iV-particle Hilbert space, that is, the subspace spanned by all 
iV-fold anti-symmetrized products of 2M single-particle orbitals (the factor 2 is due to the 
two possible spin-projections). Commonly, the many-body basis states are taken to be either 
Slater determinants (SD) or configuration state functions (CSF). The latter are special linear 
combinations of Slater determinants, which are not only eigenfunctions of the projected-spin 
operator S z - as Slater determinants are - but also eigenfunctions of the total-spin operator 
S 2 |12| . This ansatz for the wavefunction is inserted into the Schrodinger equation to obtain a 
linear equation of motion for the time-dependent expansion coefficients, which may be solved by 
a couple of methods (e.g. short iterative Lanczos- [E3H1], Chebycheff- [15], split operator- [T6] 
methods, etc. [IZ])- However, this simple approach suffers from a main drawback, namely that 
the size of the N-particle Hilbert space basis grows exponentially with the number of particles 
N and orbitals 2M, what restricts the applicability to rather small systems. A common way 
around this problem is to drop several Slater-determinants, which are believed to be physically 
less important, as it is done for instance in the time-dependent Configuration Interaction singles 
(TDCIS) method [18]. This enables the treatment of higher particle numbers resp. the inclusion 
of a larger single-particle basis at the cost of a reduced description of correlation effects. 

The multiconfigurational time-dependent Hartree-Fock method uses an alternative strategy. 
It also approximates the wavefunction by a linear expansion of basis states of H(2M, N), which 
are, however, allowed to vary in time. This is achieved by assuming the 2M orbitals { | 4>k ) } to be 
explicitly time-dependent, and being expressed by an expansion in a set of Nb time- independent 
orbitals {| \i )}'■ 

\<t>k)(t) = \ X i), k = l,-- - ,2M. (2) 
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By using Slater determinants built with such time-dependent orbitals, it is likewise only possible 
to represent states in a subspace H(2M, N). However, now this space is allowed to vary freely 
in the much larger subspace H(2Nf„N). In this way it is possible to arrive at an accurate 
description of the wavefunction and to defer the problem of the exponential growth of the 
Hilbert basis size to larger systems. 

In the following, we give the MCTDHF equations for systems described by the standard 
electronic Hamiltonian Q, which in second quantization reads 

H(t) = ^ ^ h pq (t) Epq + — y ^ Qpqrs &pqrs ; (3) 
pq pqrs 

with the one- and two-particle excitation operators |12| 

E 

J^pq 
Cpqrs 

acting in the time-dependent basis { | • The electron integrals are given by 

h pq (t) = J dr^(r)|-iA + y(r)-f(t).rU ? (r), (6) 
Qpqrs = J j drdr 0*(r)^(r) * <p*(r)(j) s (f) . (7) 

The single particle Hamiltonian h(t) includes the action of an external electromagnetic field £ (t) 
in dipole approximation (length gauge) and thus carries the only explicit time-dependence. 
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2.2. The MCTDHF equations 

As mentioned above, the multiconfigurational time-dependent Hartree-Fock ansatz approximates 
the wavefunction by a linear combination of time-dependent basis states of the iV-particle Hilbert 
space, which are in the following assumed to be Slater determinants: 



^ C n (t) | ni a ,nii3,n2a ■ ■ ■ ,riMj3;t) 



(8) 



The Slater determinants are written in occupation number representation specifying the 
occupation of the M time-dependent spatial orbitals {|0fc)} with an electron with spin- 
projection a (spin-up) or (3 (spin-down), and Ylk n ka + = N. This corresponds to a spin- 
restricted treatment, i.e. we assume that a- and /3-electrons share a common spatial orbital |19j . 
For the derivation of the equations of motion we follow Refs. [20; 21j and employ the Lagrange 
formulation of the time-dependent variational principle, in which the action functional 
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is minimized with respect to the variational parameters. The time-dependent Lagrange 
multipliers are introduced to ensure the orbitals to remain orthonormal during the temporal 
evolution. 



We first derive the equations of motion for the orbitals by requiring the variation to be 
stationary, 
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After a few steps we arrive at the following nonlinear equation: 
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where we introduced the (spin-restricted) one- and two-particle density matrices 
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as well as the mean- field operator g rs , which in coordinate representation reads 

g rs (r) = I dr' (f>*(r') - — — r s (r') . 
J \r-r'\ 

Further, the elimination of the Lagrange multipliers led to a projection operator 

P = 1 - /, | 0m)(0w | , 
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which projects on the orthogonal complement of the span of the orbitals. In order to remove the 
projection operator on the lhs. of Eq. ( 11 ) and obtain explicit equations, a unitary transformation 
among the orbitals is applied, which ensures 

0_ 
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i.e. the change of an orbital is orthogonal to the subspace spanned by all orbitals. After insertion 

(17) 



into Eq. (11), we obtain the MCTDH orbital equations 
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The minimization with respect to the coefficients, i.e. 
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then straightforwardly leads to a Schrodinger equation in matrix representation in the SD-basis: 







(19) 



Note, that again Eq. (16) has been used, which also causes the matrix element of the time- 

(20) 



derivative operator between Slater determinants to vanish 

m\ = 0. 
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Figure 1. Schematic view of the sine DVR functions, for a basis size of Nf, = 6. The 
basisfunctions are constructed over an equidistant grid in such a way, that Xi( x k) = ^ik/y/^k 
holds, with a set of integration weights W/-- Matrix elements may then be evaluated 
approximately by a summation, i.e. f dxx*( x )f( x )Xj(x) — > E/b w k X*( x k)f( x k)Xj( x k}- This 
leads directly to diagonal spatial matrix elements, ( Xi \ f( x ) \ Xj ) = f( x i) 



2.3. Numerical implementation 

We give some notes on our numerical solution of the coupled set of MCTDHF equations, Eqs. ( 19 ) 



and (|17j). 

Single-particle basis. First of all, like in Eq. ([2]), an appropriate time-independent single particle 
basis {| Xl)} 1S chosen, which inserted into the orbital equation (17) yields an equation for the 
time-dependent expansion coefficients b^^t) [19J. In this work we use a sine discrete variable 
representation (DVR) basis [22j [23] , see Fig. [I] As is common to all DVR bases, matrix elements 
of spatial operators are diagonal, and they are simply given by the function values on a grid 
of Gaussian integration points, 

(x P \f(r)\ Xq ) = S pg f(r q ), (21) 
(x P Xr\g(?,r')\xsXq) = <W $rs g(*i, rj) ■ (22) 

For the sine DVR, the grid consists of equally spaced nodes and the non-diagonal kinetic energy 
matrix can be evaluated analytically^] 

Density matrices and Hamiltonian. For the evaluation of the density matrices and the 
Hamiltonian, the matrix elements of one-and two-particle excitation operators in the basis of 
Slater determinants, ( n | E pq | m ) and ( 

| &pqrs I )> have to be evaluated, both of which may 
attain either zero, plus or minus one. Since these quantities are needed very frequently, all 
non-zero contributions are determined and their sign is stored in memory before the actual time 
propagation. The Hamiltonian can then be easily calculated using the Slater-Condon rules |12j . 
Time evolution. After writing the wavefunction coefficients C and the single-particle basis 
expansion coefficients b in a single vector, the MCTDHF equations may be casted to the general 
form 

i^v = hv), v = (£) ■ ( 23 ) 

We solve this coupled set of equations by means of general purpose integrators, such as Runge- 
Kutta or Burlisch-Stoer methods |26j. Other possible techniques particularly well suited for the 
MCTDHF scheme are given in Ref. [23j . 

Solution effort. As mentioned above, the basis of Slater-determinants grows with ( 2 ^), where 
M is the number of time-dependent spatial orbitals and the number of particles, leading 



2 Instead of a normal DVR, one could also employ a finite-element DVR, see e.g. |24l 125] . 



to the typical exponential problem of Configuration Interaction. However, for the two-particle 
model we consider here, the SD-basis only grows with 0(M 2 ) and thus poses no difficulties. 
Here, the effort is rather determined by the number jtyj, of orbitals, which may become large 
in order to adequately describe the continuum. For each evaluation of the rhs. J""(V) of 
Eq. (23), the electron integrals in the time-dependent basis are needed, which are formed 
through a transformation from the time-independent basis using the orbital coefficients. The 
time-consuming part is the transformation of the two-particle interaction matrix. Assuming 
N[> S> M, the leading term for a DVR basis is given by 0{MN?). To further diminish the effort 
and obtain an almost linear scaling with the size of the underlying time- independent basis, low- 
rank approximations of the interaction potential can be used |23j [27]. We also note, that for a 
fixed basis size Nf, the effort grows as C(M 4 A^). 



3. Numerical results 

3.1. One- dimensional Helium model 

The here considered one-dimensional model of Helium is given by a potential 

v(x) " -vTO- (24) 

and has been well tested for roughly 30 years. Compared to the real three-dimensional Helium 
atom, in this model the electron movement is restricted only to the laser polarization axis. Its 
usefulness derives on the one hand from the fact, that electronic correlation effects of Helium and 
also larger atoms can be qualitatively well explained, most prominently the nonsequential double 
ionization |28[ [29]. On the other hand, it is exactly solvable by viewing the one-dimensional 
two-particle problem as a Schrodinger equation of one-particle moving in the two-dimensional 
external potential 

2 2 1 

V(x,y) = ; + — (25) 

V^+l Vy^+l y/(x-y)* + l 

In this paper, we use it as a benchmark for the approximate MCTDHF simulations. For the 
details of the employed solution method of the two-dimensional TDSE, see [30] . 



3.2. Groundstate results 

In the following we present the results for the groundstate of the Helium model obtained with 
MCTDHF. A similar investigation has already been given in [31] . however with a slightly different 
atomic potential. The groundstate is obtained by propagation of the MCTDHF equations in 
imaginary time (ITP), starting from the orbitals of the noninteracting system. The system is 



M 


1 


2 


3 


4 


5 


exact 


Energy [Hartree] 


-2.2242 


-2.2365 


-2.2381 


-2.2382 


-2.23825 


-2.23826 


% of corr. energy 




87% 


98.9% 


99.6% 


99.9% 


100% 



Table 1. Total energies of the groundstate for different numbers M of time-dependent orbitals, 
plus the fraction of the included correlation energy (the entire correlation energy is defined as 
the difference between the exact and the Hartree- Fock result). With increasing M, the energy 
rapidly approaches the exact energy as obtained by a direct solution of the Schrodinger equation. 
For M > 5, the results equal the exact result for the given number of digits. 



Figure 2. Logarithmic contour plot of the groundstate two-particle densities for different 
numbers M of MCTDHF orbitals. x\ and X2 are the coordinates of the two electrons. In each 
plot, the red dashed curves show the result from a direct solution of the Schrodinger equation. 



evolved until the difference of the total energy between two steps falls below a certain limit 
(here lCP 10 a.u.). The results are compared to direct solutions of the Schrodinger equation, 
which were obtained by ITP as well. In Tab. [I] we display the groundstate energies for different 
MCTDHF approximations, and the corresponding fraction of the correlation energy as defined 
by the difference between the exact and the Hartree-Fock (M = 1) energy. We observe a rapid 
convergence with the number of orbitals, already the first correction to Hartree-Fock, M = 2, is 
able to account for 87% of the entire correlation energy. Note that M = 2 only needs roughly 
an eight-fold effort as compared to Hartree-Fock, which is due to the CPU time growth with 
0(M 4 N b ) (see section 

In Fig. [2]we plot the two-particle (tp) densities, that for a two-electron system equal the absolute 
square of the wavefunction, again for different approximations. In each plot, the (red) dashed 
curve depicts the tp-density from the exact solution. From the figure it is obvious, that the 
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Figure 3. Linear response of the one-dimensional Helium model atom as obtained from 
MCTDHF and the TDSE (shifted by 10 4 ). The pictures show results for different numbers 
M of time-dependent orbitals. The peaks in the tail of the spectrum are caused by correlations. 
The numbers on the upper x-axes denote the positions of the first peak, which with increasing 
M converge towards the exact position uj = 0.533 a.u. 



Hartree-Fock approximation is not able to reproduce the correct butterfly shape caused by the 
Coulomb repulsion (which is largest for x\ = x%). With increasing number of orbitals, the 
tp-density approaches the exact result, until it attains an accuracy of seven digits for M = 7. 



3.3. Time- dependent results 

As an example for time-dependent calculations, we investigate the reaction of the model atom 
to an external perturbation in linear respons^J To this end, we disturb the groundstate by 
a dipole kick, which is sufficiently short in order to provide a homogeneous spectral density 
and sufficiently weak to avoid non-linear effects. In this work, we chose a kick with duration 
0.01 a.u. and an amplitude of £q = 0.01 a.u. The disturbed system was propagated for 2000 a.u. 

3 The linear response can also be obtained from time-independent calculations by diagonalization of the 
Hamiltonian. In fact, the following results present an alternative for the diagonalization called "filter- 
diagonalization" , see e.g. [32] . 



(~ 50 fs) and the expectation value of position has subsequently been Fourier-transformed 
(using a Blackman-Harris window). This yields the dipole spectrum, from which the excitation 
frequencies can be observed. In Fig. [3] we plot the obtained spectra for different approximations 
and compare them with the result from the TDSE. First let us consider the energy region 
below the first ionization threshold uj < 0.75 a. u. [33], which corresponds to the one-electron 
excitations. These excitations are rather well reproduced already within Hartree-Fock, only the 
position of the peaks and thus the excitation energies deviate slightly. However, for Hartree- 
Fock (M = 1), there are no peaks in the tail of the curve, which consequently are caused by 
correlations. Already in the first correlation correction, M = 2, the main peaks exist, but their 
positions deviate from the TDSE result. For M = 4, the exact spectrum is well reproduced. 
The values on the upper x-axes indicate the position of the first peak, which determine the 
difference between the groundstate energy and the energy of the first excited singlet state. We 
observe a similar trend as before, namely that M = 4 is able to reproduce the exact result of 
uj = 0.533 a.u. 

4. Summary 

We have given an introduction to the time-dependent Multiconfiguration Hartree-Fock formalism 
and discussed the main ideas of our numerical implementation. In order to investigate the 
characteristics and capabilities of the method, it has been applied to a well-known one- 
dimensional model of the Helium model. We calculated the groundstate energies and two- 
particle densities, as well as the linear response of the system, each one for different MCTDHF 
approximations, i.e. for different numbers M of time-dependent spatial orbitals. All results were 
compared to direct solutions of the time-dependent Schrodinger equation. Our investigation 
showed, that MCTDHF is well suited for the time-resolved study of this two-electron system. 
We are very confident, that a comparable level of accuracy can also be obtained for larger system 
(up to, say, N = 10). 

Our future work therefore will concentrate on time-dependent electronic correlations in larger 
model atoms, which are not accessible by direct solutions of the Schrodinger equation. 
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